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^ ' Plane Couette flow perturbed by a spanwise oriented ribbon, similar to a configuration investi- 

, gated experimentally at the Centre d'Etudes de Saclay, is investigated numerically using a spectral- 

/"*\ ■ element code. 2D steady states are computed for the perturbed configuration; these differ from the 

unperturbed flows mainly by a region of counter-circulation surrounding the ribbon. The 2D steady 
flow loses stability to 3D eigenmodes at Re c = 230, /3 C = 1.3 for p = 0.086 and Re c ~ 550, /3 C ~ 1.5 for 
p = 0.043, where /3 is the spanwise wavenumber and 2p is the height of the ribbon. For p = 0.086, 
the bifurcation is determined to be subcritical by calculating the cubic term in the normal form 
equation from the timeseries of a single nonlinear simulation; steady 3D flows are found for Re as 
' low as 200. The critical eigenmode and nonlinear 3D states contain streamwise vortices localized 

near the ribbon, whose streamwise extent increases with Re. All of these results agree well with 
experimental observations. 

5h 

o ■ i. introduction 

Oh 



> 
oo 



It is well known that, of the three shear flows most commonly used to model transition to turbulence, plane Poiseuillc 
flow is linearly unstable for Re > 5772, whereas pipe Poiseuillc flow and plane Couette flow are linearly stable for all 
Reynolds numbers; see, e.g. [1]. Yet, as is also well established, in laboratory experiments, plane and pipe Poiseuillc 
flows actually undergo transition to three-dimensional turbulence for Reynolds numbers on the order of 1000. For 
plane Couette flow, the lowest Reynolds numbers at which turbulence can be produced and sustained has been shown 
to be between 300 and 400 both in numerical simulations [2,3] and in experiments [4,5]. 

The gap between steady, linearly stable flows which depend on only one spatial coordinate and three-dimensional 
turbulence can be bridged by studying perturbed versions of Couette and Poiseuillc flow. Plane Couette flow perturbed 
by a wire midway between the bounding plates and oriented in the spanwise direction has been the subject of laboratory 
experiments by Dauchot and co-workers [6-8] at CEA-Saclay. Our goal in this paper is to study numerically the flows 
£Q . and transitions in a configuration similar to that of the Saclay experiments. 

Previous studies of plane channel flows have used a variety of approaches. We briefly review these, emphasizing 
& \ computational investigations and the plane Couette case. 

One approach is to seek finite amplitude solutions at transition Reynolds numbers and to understand the dynamics of 
transition in terms of these solutions. Finite amplitude solutions for plane Couette flow have been found for Reynolds 
numbers as low as Re = 125 by numerically continuing steady states or travelling waves from other flows: the wavy 
Taylor vortices of cylindrical Taylor-Couette flow by Nagata [9,10] and Conley and Keller [11] and the wavy rolls of 
Rayleigh-Benard convection by Busse [12]. Most recently, Cherhabili and Ehrenstein [13,14] succeeded in continuing 
planc-Poiscuille-flow solutions to plane Couette flow, via an intermediate Poiseuille-Couette family of flows. They 
showed that in proceeding from Poiseuille to Couette flow, the wavespeed of the travelling waves decreases and their 
streamwise wavelength increases, as does the number of harmonics needed to capture them. When the Couette limit is 
reached, the finite amplitude solutions are highly streamwise-localized steady states. The minimum Reynolds number 
achieved in these continuations is Re = 1500. None of these steady solutions of plane Couette flow obtained so far 
are stable. 

A second, highly successful, approach has been to study the transient evolution of linearized plane Couette flow. 
Although all initial conditions must eventually decay and the most slowly decaying mode must be spanwise invariant 
by Squire's theorem, the non-normality of the evolution operator allows large transient growth. Butler and Farrell 
[15] showed that a thousand- fold growth in energy could be achieved from an initial condition resembling streamwise 
vortices which are approximately circular and streamwise invariant. Reddy and Henningson [16] computed the max- 
imum achievable growth for a large range of Reynolds numbers. An interpretation is given by these authors and by 
Trefethen et al. [17] in terms of pseudospectra: the spectra of non-normal operators display an extreme sensitivity to 
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perturbations of the operator. Thus, slightly perturbed plane Couette or Poiseuille flows may be linearly unstable for 
much lower Reynolds numbers than the unperturbed versions. 

A third broad category of computational investigation is the study of nonlinear temporal evolution in relatively 
tame turbulent plane channel flows. Orszag and Kclls [18] and Orszag and Patera [19] showed that finite amplitude 
spanwise-invariant states of plane Poiseuille flow are unstable to 3D perturbations; this is also true of quasi-equilibria 
for plane Poiseuille and Couette flow. Lundbladh and Johansson [2] showed that turbulent spots evolved from 
initial disturbances resembling streamwise vortices if the Reynolds number exceeded a critical Reynolds number 
between 350 and 375. Numerical simulations by Hamilton, Kim and Waleffe [3] of turbulent plane Couette flow 
at Re — 400 indicated that streamwise vortices and streaks played an important role in a quasi-cyclic regeneration 
process. Coughlin [20] used weak forcing to stabilize steady states containing streamwise vortices and streaks. These 
became unstable and underwent a similar regeneration cycle when the forcing or Reynolds number was increased. The 
critical Reynolds numbers displayed in all of these numerical simulations arc in good agreement with experiments by 
Tillmark and Alfrcdsson [4] and by Daviaud et al. [5] who reported turbulence at Re > 360 and Re > 370, respectively. 

The last approach we discuss, and the most relevant to this study, is perturbation of the basic shear profile, 
to elicit instabilities that are in some sense nearby. If a geometric perturbation breaks either the streamwise or 
spanwise invariance of the basic profile, then the flow is freed from the constraint of Squire's theorem, which would 
otherwise imply that the linear instability at lowest Reynolds number is to a spanwise invariant (2D) eigenmode. A 
perturbed flow with broken symmetry may directly undergo a 3D linear instability. One can hope to understand 
the behavior of the unperturbed system by considering the limit in which the perturbation goes to zero. For some 
time, experimentalists [21] have used perturbations to produce spanwise-invariant Tollmicn-Schlichting waves arising 
subcritically. More recently, for example, Schatz et al. [22] inserted a periodic array of cylinders in a plane Poiseuille 
experiment to render this bifurcation supercritical. In plane Couette flow, Dauchot and co-workers at Saclay [6-8] 
found that streamwise vortices could be induced for Reynolds numbers around 200 when a wire was placed in the flow 
(the exact range in Reynolds number for which the vortices occur depends on the radius of the wire). They suspected 
that these vortices arise from a subcritical bifurcation from the perturbed profile, but did not determine this. 

In this paper, we numerically study the destabilization of plane Couette flow when a ribbon is placed midway in the 
channel gap (Fig. 1). The ribbon is infinitely thin in the streamwise (x) direction, occupies a fraction p of the cross- 
channel (y) direction, and is infinite in the spanwise (z) direction. This geometry is similar, though not identical, to 
that used in the Saclay experiments. In the experiments, the perturbation is a thin wire with cylindrical cross-section. 
Here we use a ribbon because it is much easier to simulate numerically. For the experimental or numerical results to 
be of wider importance, the particular shape of the perturbation should not be important, as long as it is small. 

We shall address the extent to which a small geometric perturbation of the plane Couette geometry affects the 
stability of the flow. We will show that a small geometric perturbation does indeed lead to a subcritical bifurcation 
to streamwise vortices, at Reynolds numbers and wavenumbers which agree well with the Saclay experiments. 




X 

FIG. 1. Flow geometry considered in the paper. The upper and lower channel walls are separated by distance 2h and 
move with velocities Uox and — Uqx, repectively. An infinitely thin ribbon (bold line) is located midgap in the channel and has 
height 2ph (p = 0.086 for the case shown). The computational mesh (macro elements) used in our calculations is shown, as is 
the (fine) collocation mesh for polynomial order N — 8 (three elements in the enlargement). The ribbon is formed by setting 
no-slip boundary conditions on the edges of two adjoining elements. Periodic boundary conditions are imposed over length 
2Lh in the horizontal direction. The full geometry shown has aspect ratio L — 10. The system is homogeneous in the spanwise 
( z) -direct ion normal to the figure. 
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II. NUMERICAL COMPUTATIONS 



The computations consist of three parts: (1) obtaining steady 2D solutions of the Navier-Stokes equations, (2) 
determining the linear stability of these solutions to 3D perturbations, and (3) classifying the bifurcation via a 
nonlinear stability analysis. Here we outline the numerical techniques for carrying out these computations. 

A. 2D steady flows 

Our computational domain has been shown in fig. 1. We non-dimcnsionalizc lengths by the channel half-height h, 
velocities by the speed Uq of the upper channel wall, time by the convective time h/Uo. There are two nondimensional 
parameters for the flow, which we take to be the usual Reynolds number for plane Couette flow, Re = hUo/iy, where 
v is the kinematic viscosity of the fluid, and the nondimensional half-height of the ribbon p, hereafter called its radius 
for consistency with the Saclay experiments. We view the (nondimensionalized) streamwise periodicity length 2L 
as a numerical parameter which we take sufficiently large that the system behaves as though it were infinite in the 
streamwise direction. 

The fluid flow is governed by the incompressible Navier-Stokes equations: 

= _( u • V)u - Vp + — V 2 u in ft, (la) 
ut lie 

V • u = in ft, (lb) 

subject to the boundary conditions: 

u(x-L,y) = u(x + L,y) (2a) 

u(x,y = ±l) = ±x (2b) 

u(x = 0,y) = 0, for -p<y < p, (2c) 

where u = (it, v, w) is the velocity field, p is the nondimensionalized static pressure and ft is the computational 
domain. The pressure p, like u, satisfies periodic boundary conditions in x. 

Time-dependent simulations of these equations in two dimensions (w = 0, d/dz = 0) are carried out using the 
spectral element [23] program Prism [24,25]. In the spectral element method, the domain is represented by a mesh 
of macro elements as shown in Fig. 1. The channel height is spanned by five elements while the number of elements 
spanning the streamwise direction depends on its length: 24 elements are used for L — 32 and 36 elements for L — 56. 
The no-slip condition (2c) is enforced by setting zero velocity boundary conditions along the edges of two adjoining 
mesh elements: this interface defines the ribbon. If continuity were imposed along this interface, as is done on all 
other element boundaries, then the flow would reduce to unperturbed plane Couette flow. Thus the ribbon is modeled 
by a small (but significant) change in the boundary conditions on just two edges of elements in the computational 
domain. Within each element both the geometry and the solution variables (velocity and pressure) are represented 
using iVth order tensor-product polynomial expansions. The collocation mesh in Fig. 1 (enlargement) corresponds to 
an expansion with N = 8. 

A time-splitting scheme is used to integrate the underlying discrctized equations [26]. Based on simulations with 
polynomial order TV in the range 6 < N < 12 and timesteps At in the range 10~ 3 < At < 1CP 2 we have determined 
that N — 8 and At = 0.005 give valid results over the range of Re considered. These numerical parameter values 
(typical for studies of this type) have been used for most of the results reported. Each velocity component is thus 
represented by about 7500 scalars for L = 32. 

Steady flows used for our stability calculations have been obtained from simulations with Reynolds numbers in the 
range 100 < Re < 600. In all cases, the simulations were run sufficiently long to obtain asymptotic, steady velocity 
fields. We shall denote these steady 2D flows by U(x,y). 

B. Linear stability analysis 

Let \J(x,y) be the 2D base flow whose stability is sought. An infinitesimal three-dimensional perturbation 
u'(x,y, z,t) evolves according to the Navier-Stokes equations linearized about U. Because the resulting linear sys- 
tem is homogeneous in the spanwise direction z, generic perturbations can be decomposed into Fourier modes with 
spanwise wavenumbers (3: 



3 



u'(x, y, z, t) — (u cos 0z, v cos 0z, w sin 0z) 
p'(x,y,z,t) = pcos0z 



(3) 



or an equivalent form obtained by translation in z. The vector u(x,y,t) — (u,v,w) of Fourier coefficients evolves 
according to: 



= -(u • V)U - (U ■ V)u - (V - 0z)p + — (V 2 - /? 2 )u in SI, 
ut He 



(V + 0i) ■ u = 



in 17, 



(4a) 
(4b) 



where V, etc. are two-dimensional differential operators. Equations (4) are solved subject to homogeneous boundary 
conditions: 



u(x - L,y) 
u(x,y = ±1) 
u(x = 0,2/) 



u(x 


0, 



L,y) 

for -p<y<p, 



(5a) 
(5b) 
(5c) 



Equations (4) with boundary conditions (5) can be integrated numerically by the method described in section II A. 
For fixed 0, this is essentially a two-dimensional calculation [27,28]. After integrating (4)-(5) a sufficiently long time, 
only eigenmodes corresponding to leading eigenvalues remain. We use this to find the leading eigenvalues (those with 
largest real part) and corresponding eigenmodes for fixed values of Re and as follows. A Krylov space is constructed 
based on integrating (4)-(5) over K = 8 successive (dimensionless) time intervals of T = 5. More precisely, we calculate 
the fields u(t),u(t + T), . . . u(t + (K — I)T) and orthonormalize these to form a basis v\, V2, ■ ■ ■ Vk- We then define the 
K xK matrix = (vt , Lvj } where L is the operator on the right-hand-side of the linearized Navier-Stokes equations 
and (, ) is an inner product. Approximate eigenvalues a and eigenmodes u(x, y, z) are calculated by diagonalizing H; 
their accuracy is tested by computing the residual r = \ \uu — Lu||. If the eigenvalue-eigenmode pairs do not attain 
a desired accuracy (r < I0~ 5 for the case here), then another iteration is performed. The new vector is added to 
the Krylov space and the oldest vector is discarded. This is effectively subspace iteration initiated with a Krylov 
subspace. More details can be found in [22,27,29]. 

We conclude this section by considering the effect of the streamwise periodicity length 2L on the computations. 
Recall that we view L as a quasi- numerical parameter in that we seek solutions valid for large L. Figure 2 shows the 
dependence of the leading eigenvalue a on streamwise length at Re — 250, — 1.3 (values near the primary 3D linear 
instability). It can be seen that for L > 32 the eigenvalue is independent of L. This is consistent with the structure of 
the base flow and eigenmodes shown in the following section. Most of the computations reported have used L = 32. 
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FIG. 2. Leading eigenvalue as a function of streamwise periodicity half-length L for Re = 250, (5 = 1.3. For L > 32 the 
eigenvalue is independent of L. 



C. 3D simulations 



For the nonlinear stability analysis and for obtaining steady 3D flows, we carry out 3D simulations of (l)-(2) 
using the same spectral element representation in (x, y) described above together with a Fourier representation in the 
spanwise direction z. We impose periodicity in the spanwise direction by including wavenumbers m0 c for integers 
\m\ < M/2, where C is the critical wavenumber found in the linear stability analysis. The simulations we report use 
M = 16. 
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III. RESULTS 



A. 2D Steady flows 

A typical steady 2D flow for the perturbed Couette geometry is shown in Fig. 3. It is representative of base flows 
for Reynolds numbers on the order of a few hundred with a ribbon of size p = 0.086. This was chosen to correspond 
to the radius of one of the cylinders used in the Saclay experiments [6-8] . The Reynolds number Re — 250 of the 
flow shown is close to the threshold for the 3D instability that will be discussed in the next section. Unless otherwise 
stated, results are for p = 0.086, Re = 250, and L = 32. 

In Fig. 3(a) it can been seen that, except near the ribbon, the steady flow is essentially the parallel shear of 
unperturbed plane Couette flow. The streamlines are as reported experimentally in [6]. As noted there, the Reynolds 
number based on the radius of the ribbon and the local velocity near the ribbon is very small compared to Reynolds 
numbers where separation or vortex shedding could be expected. In Fig. 3(b) we plot the streamfunction of the 
deviation U — Uc where Uc = yx is the unperturbed plane Couette profile. The primary effect of the ribbon is 
to establish a region (|x| < 3) of positive circulation (opposing that of plane Couette flow) surrounding the ribbon. 
Figure 3(c) shows U — Uc over a larger streamwise extent. Further from the ribbon are wider regions (3 < \x\ < 24) 
in which the deviation is weak, but has the same negative circulation as plane Couette flow. 






FIG. 3. The steady two-dimensional base flow U(x, y) at Re — 250 for a ribbon with p = 0.086. Only the central portion of 
the full L = 32 domain is shown, (a) Streamfunction contours of U. The flow is nearly identical to the parallel shear of plane 
Couette flow except very near the ribbon, (b) Streamfunction contours of the deviation U — Uc highlighting the difference 
between the perturbed and unperturbed Couette flows. A region (|x| < 3) of positive circulation is established around the 
ribbon. The flow is centro-symmetric. (c) The deviation over a larger streamwise extent showing regions (3 < \x\ < 24) further 
from the ribbon whose circulation is negative, like that of Uc- The flow is very weak; contours of the dominant part of this 
flow are not shown. The slight lack of centro-symmetry is a graphical artifact. 

The size of the counter-rotating region is remarkably uniform over the ribbon radii p and Reynolds numbers Re 
that we have studied. We define the streamwise extent of the counter-rotating region as delimited by ip(x, y = 0) = 0, 
i.e. the x values at which the streamfunction at midhcight y = has the same value as at the channel walls y = ±1. 
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For p = 0.086, the counter-rotating region varies from |a;| < 2.18 for Re = 150 to \x\ < 3.00 for Re = 300, while 
for p = 0.043 the counter-rotating region varies from \x\ < 2.10 for Re = 150 to \x\ < 2.87 for Re = 600. This 
insensitivity to the size of p is significant in light of the 2D finite-amplitude steady states calculated by Cherhabili 
and Ehrenstein [13,14]. The states found by these authors in unperturbed plane Couette flow strongly resemble that 
in Fig. 3. These too have a central counter-rotating region surrounded by larger regions of negative circulation. At 
Re = 2200, the counter-rotating region in their flow occupies \x\ < 2.31 (see Figs. 10 and 11 of [13], Figs. 2 and 3 
of [14]) The similarity between the 2D flows for p = 0.086, p = 0.043, and, effectively, p — leads us to hypothesize 
that our 2D perturbed plane Couette flows are connected (via the limit p — > 0) to those computed by Cherhabili and 
Ehrenstein. 

We may also quantify the intensity of the counter-circulation. One measure is the maximum absolute value of v, 
which is attained very near the ribbon, at (x,y) = (±0.081,0). This value is approximately independent of Reynolds 
number, but decreases strongly with ribbon radius: v max ~ 0.031 for p = 0.086 and v max « 0.013 for p = 0.043. 

An important qualitative feature of the flow can been seen in Figs. 3(b) and (c): the flow is centro-symmetric, 
i.e. it is invariant under combined reflection in x and y, or equivalcntly rotation by angle tt about the origin. It 
can be verified that the governing equations (1) and boundary conditions (2) are preserved by the centro-symmetric 
transformation: 

u(x,y) -> -u(-x,-y). (6) 

The unperturbed plane Couette problem is also centro-symmetric. It is in fact symmetric under the Euclidean group 
Ei of translations and the "reflection" consisting of the centro-symmetric transformation (6). The ribbon in the 
perturbed flow breaks the translation symmetry, but leaves the centro-symmetry intact. Note that reflections in x or 
y alone are not symmetries of either the unperturbed or the perturbed plane Couette problem because either reflection 
alone reverses the direction of the channel walls, violating the boundary conditions (2b). 

In Fig. 4 we present streamwise velocity profiles near the ribbon. For \x\ > 0.5, the Couette profile is very nearly 
recovered. Figure 4(b) shows streamwise velocity profiles of the deviation from the linear Couette profile across the 
full channel. Close examination reveals that these profiles are not odd in y, consistent with the fact that the system 
is neither symmetric nor antisymmetric under reflection in y. The symmetric partners to the profiles shown are at 
negative x values. 

The profiles in Fig. 4 are similar to those Bottin et al. [8] obtained in the Saclay experiments under similar 
conditions. It is not possible to compare directly with experiment because of the difficulty in obtaining experimental 
velocity profiles and because the geometric perturbations differ in the computations and experiments. The only 
noticeable difference between experiments and computations is that the profiles Fig. 4(b) are very nearly odd in y, 
whereas in experiment this lack of symmetry is more pronounced. 




FIG. 4. Streamwise velocity profiles in the perturbed geometry, (a) U(x,y) as a function of y for x = 0.081 (dotted), 
x = 0.25 (dashed), and x = 0.5 (solid). Only the central portion of the channel is shown. Only very close to the ribbon does 
the velocity differ significantly from the linear profile, (b) Deviation U(x,y) — y over the full range of y. 
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Finally in Fig. 5 we quantify the deviation between perturbed and unperturbed plane Couette flow by plotting the 
energy per unit length E(TJ — Uc) = j\ dy^\U(x,y) — Uc(j/)| 2 as a function of x for —56 < x < 56. The data show 
a narrow central region, corresponding to the region \x\ < 2.76 of positive circulation seen in Fig. 3(b), where the 
deviation falls sharply and approximately exponentially in x. For \x\ > 2.76, the deviation, while very small, decays 
very slowly (and not exponentially) with \x\. The boundaries x = ±23.78 terminating the outer region of negative 
circulation can also be seen on fig. 5. The precision of the computations is surpassed beyond \x\ = 40. This figure 
shows that for \x\ > 32 the deviation of the base flow from Couette is indeed very weak and this supports our choice 
of L = 32 as an adequate domain size for most computations. 
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FIG. 5. Energy of deviation between perturbed and unperturbed Couette flows as a function of x. Parameters are the same 
as in Fig. 3 except that here L — 56. Abrupt changes in slope at \x\ — 2.76, \x\ = 23.78 correspond to changes in the sign of 
the circulation of U — Uc- For |a;| > 40, the deviation is below the precision of the computations. 

B. 3D Linear stability results 

The two-dimensional steady flows just discussed become linearly unstable to three-dimensional perturbations when 
the Reynolds number exceeds a critical value Re c . To determine this value and the associated wavenumber, we have 
performed a linear stability analysis of the steady flows via the procedure described in Sec. II B. 

Figure 6 shows the growth rate a of the most unstable three-dimensional eigenmode u as a function of Re and 
spanwise wavenumber /3 for a ribbon with p = 0.086. For each value of Re, we have fit a piecewise-cubic curve, 
shown in fig. 6, through the eigenvalue data to determine the wavenumber (3 max (Re) which maximizes a. The critical 
Reynolds number Re c is then determined by linear interpolation of cr(/3 max (i?e), Re) through these maxima and finding 
its zero crossing. From this we find critical values for the onset of linear instability to be Re c = 230 and fi c = 1.3 
for the ribbon with p = 0.086. These values are consistent with what is seen experimentally, but we delay discussion 
until the conclusion. 

Figure 7 shows similar eigenvalue spectra for a ribbon half as large: p — 0.043. The critical wavenumber fi c w 1.5 
is only slightly larger than the previous value. However, the critical Reynolds number is much larger: Re c w 550. 
The critical Reynolds number must increase as p is decreased since, when no ribbon is present, the problem reduces 
to classical plane Couette flow which is linearly stable for all finite Re, i.e. lim p ^o Rgc(p) = oo. 

We note that Cherhabili and Ehrenstein [14] also calculate 3D instability for their 2D finite amplitude plane Couette 
flows. Despite the resemblance of their 2D flows to ours, the spanwise wavenumber corresponding to maximal growth 
is much larger in their case: f3 k, 23. 
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FIG. 6. Growth rate a of most unstable three-dimensional eigenmode as a function of spanwise wavenumber f3 for 
Re = 150, 200, 250, 300 with ribbon radius p = 0.086. Critical values for instability are Re c = 230 and f3 c = 1.3. 
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FIG. 7. Growth rate of most unstable eigenmode as a function of j3 and Re for p = 0.043. Critical values are Re c « 550 and 
/3 C « 1.5. 

A computed eigenvector u = (u, v, w) is shown in Figs. 8 and 9. This eigenvector is near-marginal: Re — 250, close 
to Re c — 230. The spanwise wavelength is A = A c = 2ir/ j3 c = 4.83. The other parameters are p — 0.086 and L = 32. 

Figure 8 shows (v, w) velocity plots at four streamwise locations. In the x = plane containing the ribbon, the flow 
is reflection-symmetric in y, and the flow is primarily spanwise. The trigonometric dependence in z with the choice 
of phase (3) can be seen. In the planes x = 1, x = 2, and x — 3, two counter-rotating streamwise vortices are present. 
The flow for negative x is obtained by reflection in y. 
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FIG. 8. Velocity field of near- marginal eigenvector in the planes x = 0, x = 1, x = 2, and x = 3. Parameters are 7?e = 250, 
p = 0.086, L = 32, A = 4.83. At x — 0, the velocity is reflection-symmetric in y and primarily spanwise. The ribbon is seen as 
the area \y\ < p = 0.086 with no flow. Streamwise vortices are visible for x > 1. The scale for distances in x is stretched by a 
factor of 3.333 relative to distances in y,z. 




FIG. 9. Velocity field of near-marginal eigenvector in the planes z = and y = 0. 

Figure 9 presents two complementary views of the eigenvector ufor— 1 < x < 13. Above is a plot of (u,v) in the 
plane z = where they are maximal (cf. eq. 3). Below is a plot of (u, to) in the plane y — at mid-channel height. 
Figure 10 shows the x— dependence of the spanwise-averaged energy per unit length 
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E{u) = ±-£dz y^dy ||u| 2 (7) 

Here, the eigenvector was computed in a larger domain (L = 56) in order to determine its long-range behavior. The 
eigenvector is localized: the energy decays exponentially with |x| and does not reflect the counter- and co- rotating 
regions of the 2D base flow seen on Fig. 5. The flow deficit due to the ribbon produces the local minimum at x = 0. 
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FIG. 10. Energy of near-marginal eigenvector as a function of x. Parameters are the same as in Figs. 8 and 9 except that 
here L — 56. Vertical scale is arbitrary. 

We have also computed the vorticity of the eigenvector. Despite the streamwise vortices visible on Fig. 8, u x is the 
smallest vorticity component and u z by far the largest over most of the domain. 



C. 3D Nonlinear stability results 

Our method of nonlinear stability analysis has previously been used to determine the nature of the bifurcation to 
three-dimensionality in the cylinder wake [28]. The method is based on tracking the nonlinear evolution of the 3D 
flow starting from an initial condition near the bifurcation at Re c . "Near" refers both to phase space (i.e. a small 3D 
perturbation from the two-dimensional profile) and to parameter space (i.e. at a Reynolds number slightly above the 
linear instability threshold). In essence we follow the dynamics along the unstable manifold of the 2D steady flow 
far enough to determine how the nonlinear behavior deviates from linear evolution. From this we can determine very 
simply whether the instability is subcritical or supercritical. 

Three-dimensional simulations are carried out for p = 0.086 at Re = 250, slightly above Re c = 230, starting with 
an initial condition of the form: 



u(x, y, z) = U(x, y) + eu(x, y, z), 



(8) 



where U is the 2D base flow at Re = 250, u is its eigenmodc at wavenumber (3 C — 1.3, and e is a small number 
controlling the size of the initial perturbation. 

The restriction to wavenumbers which are multiples of [3 C accurately captures the evolution from initial condition 
(8), since the Navier-Stokes equations preserve this subspace of 3D solutions. That is, we seek only to follow the 
evolution in the invariant subspace containing the critical eigenmode. We do not address the issue of whether the 
A c -periodic flow is itself unstable to long- wavelength perturbations. 

To analyze the nonlinear evolution, we define the (real) amplitude A of the 3D flow as: 



[A 



^ r^c /■+! r+L ^ 



1/2 



(9) 



where Ui(x, y, z, t) is the component of the 3D velocity field at wavenumber (3 C , i.e. A is the square root of the energy 
of the flow at wavenumber (3 C . (A complex amplitude, not required here, would include the phase of the solution in 
the spanwise direction.) 
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Figure 11 shows the time evolution of A from our simulations. The value of e is such that the initial energy of the 
3D perturbation, eu, is E = A 2 = 1.66 x 10~ 5 ; the energy of the base flow U is E = 21.3. 
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FIG. 11. Nonlinear growth of the amplitude A of the 3D flow from simulation (solid) at Re — 250 plotted on linear and 
logarithmic scales. First-order (dotted) and third-order (dashed) dynamics are shown with a — 0.004669 and a — 0.9. The 
faster than exponential nonlinear growth (i.e. positive a) shows that the instability at Re c is subcritical. 

To interpret the nonlinear evolution, consider the normal form for a pitchfork bifurcation including terms up to 
third-order in the amplitude: 

A = a A + aA 3 (10) 

The leading nonlinear term is cubic because the 3D bifurcation is of pitchfork type (an 0(2) symmetric pitchfork 
bifurcation). The Landau coefficient a determines the nonlinear character of the bifurcation. If a > 0, then the 
nonlinearity is destabilizing at lowest order and the bifurcation is subcritical; if a < 0, then the cubic term saturates 
the instability and the bifurcation is supercritical. 

Figure 11 includes curves for first-order evolution (i.e. A = aA) and the third-order evolution given by Eq. (10). 
For the first-order evolution, the eigenvalue a for the bifurcation has been computed via the linear stability analysis in 
Sec. IIIB. For the third-order evolution we have simply fit the one remaining parameter, a, in the normal form. We 
followed the approach in [28] of using the time series A{t) and the known value of a to estimate a from a ~ {A— a A) /A 3 . 
This gives a = 0.9 ±0.05, a constant value for T < 500, which determines how long the third-order truncation is valid 
in this case. The value of a is essentially unchanged when the mesh is refined by increasing the polynomial order iV 
to 10 or the number M of Fourier modes to 32. The magnitude of a depends on the definition of A, but its sign does 
not. The fact that a is positive indicates that the instability is subcritical. Figure 11 indicates that the 3D flow has 
become steady by t « 1000. The nonlinear saturation seen in the time series is not captured by including a fifth-order 
term in the normal form. 

We have verified that the instability is subcritical by computing nonlinear states below Re c . In Fig. 12, we show 
the steady 3D flow at Re — 200. This figure is analogous to fig. 8 depicting the eigenvector, so we will emphasize 
here the ways in which the two flows differ. Small streamwise vortices can be seen in each of the four corners of the 
x = plane containing the ribbon. The lower (y < 0) pair evolve with x into the strong pair of vortices at x = 1. The 
vortices at x — 3 are tilted with respect to their counterparts in the eigenvector, attesting to the nonlinear generation 
of the second spanwise harmonic 2/3. The 3D flow in the y = and z = planes (after subtraction of the dominant 
2D base flow) is sufficiently similar to the eigenvector (fig. 9) that we do not present it here. 
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FIG. 12. 3D velocity field in the planes x = 0, x = 1, x = 2, and x = 3. Parameters are Be = 200, p = 0.086, L = 32, A = 4.83. 
At a; = 0, four small streamwise vortices can be seen in the corners of the domain. The lower (y < 0) vortex pair evolves into 
the large vortices seen at x = 1. The scale for distances in x is stretched by a factor of 3.333 relative to distances in y,z. 




FIG. 13. Streamwise velocity contours in the plane x = 2 for the 3D field. Solid contours correspond to u > 0, dashed 
contours to u < 0. u is small and nearly constant in the interior near z = 0, A; elsewhere it varies approximately linearly with 

y- 



Streamwise velocity w-contours of the 3D flow at x — 2 are shown in Fig. 13. The (v, w) projections of our 3D 
flow in Fig. 12, showing the tilted streamwise vortices, resemble the depictions of optimally growing modes by Butler 
and Farrell [15], of instantaneous turbulent flows by Hamilton et al. [3] and of weakly forced states by Coughlin [20]. 
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However, our streamwise velocity u pictured in Fig. 13 differs significantly from [3,20] in that their u contours are 
much more strongly displaced at the vortex boundaries. This is probably due to the fact that our Reynolds number 
of 200 is substantially lower than their Re — 400. 

Far from the ribbon, the 3D flow returns to plane Couette flow. Fig. 14 compares the energy distribution E(U — Uc) 
defined by (7) of the deviation of the 3D flow from plane Couette flow at Re = 200 with that at Re = 250. Note that 
the 3D flow is less localized than the corresponding eigenvector (Fig. 10) . It can be seen that at the higher Reynolds 
number the deviation has higher energy, and importantly, occupies a larger streamwise extent. This is in accord with 
the experimental observation that the streamwise extent of vortices in the perturbed flow increases with increasing 
Reynolds number. 




FIG. 14. Energy of deviation from plane Couette flow of 3D velocity fields at Re = 200 (dashed) and Re — 250 (solid) as a 
function of x. Plotted is the square root of energy per unit length in x. The range in x is taken larger than the computational 
domain (L — 32) to match the range of Figs. 5 and 10. 

We have attempted to determine the location of the saddle-node bifurcation marking the lower Reynolds number 
limit of this branch of steady 3D states; we believe that it occurs just below Re = 200. There remains nevertheless a 
slight uncertainty regarding the lower bound for these states because we have found evidence of two different types of 
branches of steady 3D states over the range 200 < Re < 250. The study of these states is further complicated by the 
fact that the time evolution to many of them is oscillatory, indicating that their least stable eigenvalues are a complex 
conjugate pair. Further investigation is required to ascertain the full nonlinear bifurcation diagram. 

We have also sought to determine how the scenario changes as the ribbon radius p is decreased. Recall from 
section IIIB that for p = 0.043, we found Re c w 550. At these parameter values, 3D simulations display chaotic 
time evolution. By decreasing Re, we have succeeded in computing a stable 3D steady state at Re = 350. Since the 
simulation showed chaotic oscillation for a long time (3000 time units) before showing signs of approaching a steady 
state, there remains the possibility that stable 3D steady states are also attainable for higher Re. Simulations at 
Re = 300 result in decay to the basic 2D state (although we do not exclude the possibility of maintaining 3D states 
by a more gradual decrease in Re) . This is consistent with simulations of the unperturbed plane Couette geometry 
(p = 0) by Hamilton et al. [3], who observed chaotic oscillation for Re > 400 and plane Couette flow Uc = yx for 
Re = 300. Other numerical [2] and experimental [4,5] investigations in the unperturbed plane Couette geometry also 
indicate a critical Reynolds number of 360-375 for transition to turbulence. We plan to investigate the p-dependence 
of the steady 3D states and their stability in a future publication. 



IV. CONCLUSION 



We have performed a computational linear and nonlinear stability analysis of perturbed plane Couette flow in order 
to understand experiments recently performed at Saclay [6-8], and more generally, three-dimensional flows in the 
plane Couette system. We have accurately determined the extent to which the basic steady 2D profile is modified by 
the presence of a small spanwisc-oricntcd ribbon in the flow. We have determined that such a ribbon, comparable in 
size to the cylinders used in the Saclay experiments, is large enough to induce linear instability of the basic profile at 
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Reynolds numbers of order a few hundred. 

We elaborate further on how our analysis complements the Saclay results. An experimental diagram was obtained 
[7,8] for the Reynolds number range of existence of various types of flows: 2D, 3D with streamwise vortices, inter- 
mittent, and turbulent. In these experiments it was not determined whether the 3D streamwise vortices arise from a 
linear instability of the 2D flow. Our results show that a small geometric perturbation does destabilize the 2D flow 
in a subcritical instability and that the bifurcating solution is a 3D flow with streamwise vortices. Specifically, for a 
nondimensional radius p = 0.086, we find Re c = 230 and for p = 0.043, we find Re c = 550. The computed spanwise 
wavelength of the most unstable mode is in good agreement with the value seen experimentally. The streamwise 
extent occupied by these vortices decreases with decreasing Reynolds number, as observed in experiment, and is finite 
at the lower Reynolds limit of the 3D flows. 

The Reynolds number ranges for the steady 3D flows we have computed differ somewhat from those seen exper- 
imentally. For p = 0.086, streamwise vortices were observed experimentally over the range 150 < Re < 290. For 
p = 0.086, we have thus far found steady 3D flows only if Re > 200. In experiments with p = 0.043, streamwise 
vortices have been observed over the range 190 < Re < 310; we have thus far found steady 3D flows only for Re near 
350. However, a full study of the 3D flows is still pending and may resolve these discrepancies. 

There are also minor qualitative differences between our results and the experimental findings. The first is that 2D 
flows in our computations are more antisymmetric in the cross-channel direction y than in experiment (our Fig. 4 vs. 
Fig. 4 of [8]). This is probably due to the fact that we perturb our flow with an infinitely thin ribbon and not a wire 
(cylinder) as in the experiment. However, based on the existence of instability in both cases, this small difference in 
the basic 2D flow is probably not very significant. The other difference is due to the fact that our 3D simulations 
impose spanwise periodicity with a single critical wavelength A c . Experimentally, it is observed that the streamwise 
vortices are not always regularly spaced in the spanwise direction. 

As stated in the introduction, streamwise vortices can be made to appear in channel flows via a number of ap- 
proaches. Butler and Farrell [15] and Rcddy and Hcnningson [16] show that under linear evolution modes of this type 
can achieve very high amplitude before eventually being damped. In the pseudospectrum interpretation of Trefethen 
ct al. [17,16], non-normality leads to sensitivity of the spectrum: streamwise vortices are unstable modes of a slightly 
perturbed linear stability matrix. Our results are entirely consistent with this interpretation: the ribbon or wire serves 
as a specific realization of a perturbation to the stability matrix, and has indeed rendered the flow linearly unstable 
to streamwise vortices. 

Future computational work is needed to explore these flows. Calculating complete bifurcation diagrams for the two 
cases p = 0.086 and p — 0.043 is the first priority. We plan to study quantitatively and qualitatively the bifurcations 
at which the steady 3D flows terminate at low Re and lose stability at high Re. Our goal is to continue 2D and 3D 
solutions of the perturbed system to the plane Couette case, p = 0. 
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